options(scipen = 999)

# Load the necessary libraries
library(Synth)
library(dplyr)
library(gsynth)
library(writexl)

data <- read.csv("base_4dic.csv", row.names = 1) #GDP Per Capita Constant 2015 - WORLD BANK
data <- data %>% subset(year <= 2009)

data$country_name <- as.character(data$country_name)
data$country_code <- as.character(data$country_code)

# country_left = 'Todos los paises'
country_left = c('China', 'Philippines', 'Canada', 'United States', 'South Africa','Australia') # BENCHMARK

  treated_country = 'Chile'
  intervention_year = 1990 #1990(correct) - 1995 Time-Placebo - 1983 No effect on growth (Chicago boys etc...)
  
  # Identify the pre-treatment and post-treatment periods
  pre_treatment_period <- c(1975:(intervention_year-1))   # Replace with the years of your pre-treatment period
  post_treatment_period <- c(intervention_year: 2009)  # Replace with the years of your post-treatment period
  
  # Create the outcome and covariates matrices
  outcome <- data[data$country_name == treated_country, "gdp_percapita"]
  
  covariate_names <- c(names(data)[4:16]) #[5:16]
  
  # Chequeando variables excluidas
  names(data)[!names(data) %in% covariate_names]
  
  covariates <- data %>% subset(country_name != treated_country) %>% select(all_of(covariate_names))
  
  data$country_index <- as.factor(data$country_code)
  data$country_index <- as.numeric(data$country_index)
  
  # Guarda los resultados
  
  library(car)
  model <- lm(gdp_percapita ~ ., data=data %>% select(c(all_of(c('gdp_percapita', covariate_names)))))
  vif_values <- vif(model)
  print(vif_values)
  
  treated_country_index = unique(data$country_index[which(data$country_name == treated_country)]) #dejar afuera tratado
  left_country_index = unique(data$country_index[which(data$country_name == country_left)]) #dejar afuera leave one out
  if(length(treated_country_index) > 1){treated_country_index = treated_country_index[1]}
  
  input_dataprep_para_SC <-				
    dataprep(data,				
             predictors = covariate_names,				
             dependent     = "gdp_percapita",				
             unit.variable = "country_index",				
             time.variable = "year",				
             unit.names.variable = "country_name",				
             treatment.identifier  = treated_country_index,				
             controls.identifier   = c(1:24)[!c(1:24) %in% c(treated_country_index, left_country_index)], #all the 98 regions from the donor pool				
             time.predictors.prior = pre_treatment_period,				
             time.optimize.ssr     = pre_treatment_period,				
             time.plot             = c(pre_treatment_period,post_treatment_period))				
  
  
  modelo.SC <- synth(input_dataprep_para_SC)	
  
  tablas_resultados.SC   <- synth.tab(				
    dataprep.res = input_dataprep_para_SC,				
    synth.res    = modelo.SC)				
  
  # Pesos países
  w.at11 <- as.data.frame(tablas_resultados.SC$tab.w)				
  w.at11				
  
  w_suficiente <- w.at11 %>% subset(w.weights > 0.001)
  
  # Importancia relativa variables
  v.at11 <- as.data.frame(tablas_resultados.SC$tab.v)				
  v.at11				
  
  v_suficiente <- v.at11 %>% subset(v.weights > 0.001)
  
  
  #-----Matching treated vs Synth--------
  matching <- as.data.frame(tablas_resultados.SC$tab.pred)	
  matching
  
  synth.out.at11 <- synth(input_dataprep_para_SC)	
  
  synth.tables.at11   <- synth.tab(				
    dataprep.res = input_dataprep_para_SC,				
    synth.res    = synth.out.at11)	
  
  # Crear un gráfico para visualizar los resultados --------
  path.plot(
    dataprep.res = input_dataprep_para_SC, # los datos preparados
    synth.res = synth.out.at11,       # los resultados de la función synth()
    Main = "Chile's GDP Per Capita",    # el título principal del gráfico
    Ylab = "GDP per capita",         # la etiqueta del eje y
    Xlab = "Year",                    # la etiqueta del eje x
    Legend = c("Chile", "Synthetic Chile"), # la leyenda
    # CI = TRUE                        # si quieres incluir un intervalo de confianza
  )
  
  
  # Agrega una línea vertical en el año XXX (color rojo)
  abline(v = 1990, col = "red") #1990
  #abline(v = 1983, col = "red") #Placebo test - No effect, same divergence
  
  
  # Generando vector con valores de Chile Sintético--------------------------------
  
  data_synth_chile <- input_dataprep_para_SC$Y0plot%*%synth.out.at11$solution.w
  
  data_synth_chile <- cbind.data.frame(country_name = 'Synthetic Chile',
                                       year = row.names(data_synth_chile),
                                       gdp = data_synth_chile)
  
  data_synth_chile <- data_synth_chile %>% rename(gdp = w.weight)
  data_synth_chile
  
  
# Plot Synth Control -------------------------------------------------------
  
  # Convertir las columnas a los tipos de datos correctos
  chile_sintetico <- (input_dataprep_para_SC$Y0plot %*% modelo.SC$solution.w)
  
  chile_sintetico <- cbind.data.frame(year = row.names(chile_sintetico) ,
                                      gdp_percapita = chile_sintetico)
  row.names(chile_sintetico) = NULL
  names(chile_sintetico) <- c('year', 'gdp_percapita')
  
  data_chile_org_y_synth <- rbind.data.frame(cbind(chile_sintetico, pais = 'Chile sintético'),
                                             cbind(datos_chile_original, pais = 'Chile original'))
  data_chile_org_y_synth$year <- as.numeric(as.character(data_chile_org_y_synth$year))
  
  # data_total <- data_total %>% subset(country_name %in% c('Chile', 'Synthetic Chile', as.character(paises_peso_positivo$unit.names)))
  # data_chile_org_y_synth <- data_chile_org_y_synth %>% subset(pais %in% c('Chile', 'Synthetic Chile'))#, as.character(paises_peso_positivo$unit.names)))
  data_chile_org_y_synth$gdp_percapita <- as.numeric(data_chile_org_y_synth$gdp_percapita)
  
  # Generar el gráfico
  ggplot(data_chile_org_y_synth, aes(x = year, y = gdp_percapita, color = pais)) +
    geom_line(aes(alpha = 1)) +
    # geom_line(aes(alpha = ifelse(pais %in% c("Chile", "Synthetic Chile"), 1, 0.3))) +
    scale_alpha(guide = 'none') + 
    labs(x = "Year", y = "GDP per capita", color = "Legend") +
    theme_minimal() +
    geom_vline(xintercept = 1990, linetype = 'dashed', colour = 'gray')
  
  
